arXiv: 1503.0267lv5 [cond-mat.soft] 8 Oct 2015 


The Raspberry Model for Hydrodynamic Interactions Revisited. I. 
Periodic Arrays of Spheres and Dumbbells 


Lukas P. Fischer/ Toni Peter/ Christian Holm/ and Joost de Graaf^’|^ 

^Institute for Computational Physics (ICP), University of Stuttgart, Allmandring 3, 70569 Stuttgart, Germany 

(Dated: October 12, 2015) 

The so-called ‘raspberry’ model refers to the hybrid lattice-Boltzmann and Langevin molecular 
dynamics scheme for simulating the dynamics of suspensions of colloidal particles, originally devel¬ 
oped by [V. Lobaskin and B. Diinweg, New J. Phys. 6, 54 (2004)], wherein discrete surface points 
are used to achieve fluid-particle coupling. This technique has been used in many simulation studies 
on the behavior of colloids. However, there are fundamental questions with regards to the use of 
this model. In this paper, we examine the accuracy with which the raspberry method is able to 
reproduce Stokes-level hydrodynamic interactions when compared to analytic expressions for solid 
spheres in simple-cubic crystals. To this end, we consider the quality of numerical experiments that 
are traditionally used to establish these properties and we discuss their shortcomings. We show 
that there is a discrepancy between the translational and rotational mobility reproduced by the 
simple raspberry model and present a way to numerically remedy this problem by adding internal 
coupling points. Finally, we examine a non-convex shape, namely a colloidal dumbbell, and show 
that the filled raspberry model replicates the desired hydrodynamic behavior in bulk for this more 
complicated shape. Our investigation is continued in [J. de Graaf, et at, J. Chem. Phys. 143, 

084108 (2015)], wherein we consider the raspberry model in the confining geometry of two parallel 
plates. 


I. INTRODUCTION 

The physical description of hydrodynamic interactions 
in fluids has been a field of intensive study for over three 
centuries. The first mathematical description of (rarified) 
flow dates back to Euler. [1] This description was subse¬ 
quently refined by Navier and Stokes to be applicable to 
the flow of dense media. Hi However, finding solutions 
to the Navier-Stokes equations, even under the simplify¬ 
ing assumption of the low Reynolds number regime, has 
proven to be a particularly challenging boundary-value 
problem. Only in a few simple geometries can the Navier- 
Stokes equations be analytically solved, often leading to 
truncated series expansions rather than a full solution. 

Two geometries that can be handled semi-analytically 
are a simple-cubic array of spheres and a sphere between 
two parallel plates. The former is of particular interest 
as a toy model for fluid flow in a porous medium (at 
small sphere separations), [1] while the latter is relevant, 
for example, to the field of hydrodynamic chromatog¬ 
raphy. Hi In this paper, we consider the crystalline 
arrangement and in Part H of our investigation, [7] we 
study the confining geometry of two parallel plates. 

There are a myriad of (semi-) analytic investigations 
for the simple-cubic geometry, which makes this geome¬ 
try perfectly suited for benchmarking the quality of hy¬ 
drodynamic solvers. For the translational movement of a 
simple-cubic crystal through a fluid, the first results were 
obtained by Hasimoto, who derived a semi-numerical re¬ 
sult for dilute systems. [5] A complete numerical study 
for a larger range of lattice spacings and various crystal 
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structures was later presented by Zick and Homsy. [9] The 
hydrodynamic flow around an infinite (simple-cubic) ar¬ 
ray of rotating spheres was first described by Brenner et 
al. [To] These results were subsequently refined by Zu- 
zovsky et al. El A complete numerical study of both 
translational and rotational friction over a large range 
of possible lattice spacings was provided by Hofman et 
al. |4] We utilize this large body of data as a reference 
throughout our manuscript. 

A breakthrough in the numerical simulation of fluid 
dynamics resulted from the development of the lattice- 
Boltzmann (LB) algorithm. LB is based on the dis¬ 
cretized version of the Boltzmann transport equation, 
see, e.g., Ref. [T2] for a brief background. This lattice- 
based algorithm allows for the efficient simulation of hy¬ 
drodynamic interactions in arbitrary geometries using 
simple boundary conditions, such as the bounce-back rule 
to obtain no-slip surfaces. [12] 

One method to model particles moving in an LB fluid 
was introduced by Ahlrichs and Diinweg, who simu¬ 
lated polymer chains by utilizing an interpolated point¬ 
coupling scheme. [I3] These points couple to the fluid 
through a frictional force, acting both on the solvent and 
on the solute, which depends on the relative velocity. The 
effect of this coupling is the formation of a hydrodynamic 
hull around the points, which thus gain a finite hydrody¬ 
namic extent (effective hydrodynamic radius). [13] Even 
if individual friction coefficients, and thus different effec¬ 
tive radii, are used for the points, this method is limited 
in the effective size ratios that it can handle. Namely, by 
the particle-grid interpolation scheme and discretization 
used for the LB fluid. [T3] Thus the method cannot be 
employed to study systems with substantial variation in 
particle size, for example, the electrophoresis of colloids 
with explicit ions. 
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FIG. 1. (color online) Representation of the structure of the 
raspberry models used in our simulations, filled (left) and 
hollow (right), respectively. The central bead to which all 
other beads are connected via rigid bonds, is shown using a 
green sphere. The blue spheres represent the beads that form 
the filled raspberry and the red ones give the surface beads 
used for the hollow variant. The radius of the beads is chosen 
to be smaller than the typical effective hydrodynamic radius 
to help visualize the internal structure. This figure is also 
included in Part II [7] of our analysis of the raspberry model. 


Lobaskin and Diinweg remedied this issue by intro¬ 
ducing the so-called ‘raspberry’ model, in which a larger 
colloid is modeled using the aforementioned coupling by 
discretizing the surface of the colloid into points. [14] The 
method derives its name from this discretized nature of 
the surface, which resembles a raspberry, when repre¬ 
sented by molecular-dynamics (MD) beads, see Fig.[^ A 
proper coverage of the surface by coupling points, such 
that the fluid inside of the shell is ‘trapped’ and thus 
translates and rotates in unison with the shell, was as¬ 
sumed to create an effective no-slip/co-moving boundary 
condition at the surface. [HdS] 

Other methods to simulate moving boundaries exist, 
both LB-based and non-LB. Moving Ladd bounce-back 
(Ladd BB) boundaries exploit the lattice structure of 
the LB in describing colloidal particles. m dH The 
immersed boundary method (IBM), [TH] and the ex¬ 
ternal boundary force (EBF) method [IH] both use a 
point-coupling strategy to describe particles in LB; al¬ 
though it has recently been shown that these methods 
are special choices of the friction and mass ratio in the 
Ahlrichs and Diinweg scheme. m The most commonly 
used non-LB methods include: dissipative particle dy¬ 
namics (DPD), [^[55] multi-particle collision dynamics 
(MPCD) or stochastic rotation dynamics (SRD), [231 [21] 
and Stokesian Dynamics (SD). [25] However, the rasp¬ 
berry method has remained popular, because of its sim¬ 
plicity as a straightforward extension of point-particle 
coupling. It has been extended upon [IS] and has been 


used in a wide variety of simulation settings. [2M^ Re¬ 
cently, this model was employed in the context of multi¬ 
particle collision dynamics (MPCD), [531 HO] stochastic 
rotation dynamics (SRD), [31] and dissipative particle 
dynamics (DPD) simulations. [351133] 

Of singular interest are a set of recent publications 
from the Denniston group. [33H37] In these publica¬ 
tions the quality of the (raspberry-type) point-coupling 
schemes are investigated and compared to theoretical ex¬ 
pressions. Ollila et al. show in Ref. [34] that there is good 
correspondence between the LB simulations and analytic 
results [3H1I33] for a hollow shell, an annulus, and a dense 
distribution of coupling points. They place these results 
in the context of the simulation of porous particles. In 
Ref. [3S] , Ollila et al. further analyzed the quality of the 
point-coupling method and showed that there are prob¬ 
lems with this scheme when utilizing it to describe solid 
particles. In particular, Ollila et al. demonstrated that 
the hydrodynamic radius of these particles is ill-defined 
in an LB fluid. That is, the effective hydrodynamic ra¬ 
dius that follows from the translational mobility {via the 
Stokes relation) does not match that obtained using the 
rotational mobility. By careful calibration, [S^ the use 
of a colloid radius that is ‘incommensurate’ with the lat¬ 
tice spacing, [34] and modification of the coupling of the 
points to the LB fluid, [31 [37] the rotational and trans¬ 
lational effective radii can be well-matched in the for¬ 
malisms suggested by Ollila et al. and coincide with the 
radius given to the coupling points. 

In this manuscript, we re-examine the raspberry model 
by Lobaskin and Diinweg m in the context of the work 
of Ollila et al. [31 [3S] We show that there is a simple 
way to obtain an effectively consistent hydrodynamic de¬ 
scription of a solid particle using the raspberry model for 
suitably chosen LB and coupling parameters. Namely, 
by a ‘filling -|- fitting’ strategy, which we will describe 
in detail in Section (IIIA2). This approach consists of 
introducing coupling points to the interior of the rasp¬ 
berry particle and fitting for the radius of a solid par¬ 
ticle using suitable experiments. Our ‘filling -|- fitting’ 
procedure does not necessitate a particle radius that is 
incommensurate with the lattice. Moreover, it yields an 
internally consistent formalism, which reproduces the hy¬ 
drodynamic properties of a solid object with a high de¬ 
gree of accuracy. 

We show how our fit parameter (the effective hydrody¬ 
namic radius) can be straightforwardly determined. To 
demonstrate that our method works for a range of reason¬ 
able LB parameters, we examine the quality of the rasp¬ 
berry model in the classic fluid-dynamics geometry of a 
simple-cubic arrangement. diidniiii] We show that the 
raspberry model reproduces the theoretical result sur¬ 
prisingly well over the complete range of applicable rasp¬ 
berry (sphere) separations. In obtaining these results, 
we also analyze the quality of the standard hydrody¬ 
namics experiments performed in this geometry. |14[ 115] 
We further demonstrate that the improved correspon¬ 
dence between the effective rotational and translational 
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hydrodynamic radius is upheld over a large range in bare 
frictions, (original/imposed) radii of raspberry, and suf¬ 
ficiently large filling fractions. We also comment on the 
interpretation of our data in the context of theoretical 
results for porous objects. Finally, we consider 

the effectiveness of the raspberry description in model¬ 
ing solid non-convex particles and show that the ‘filling 
-|- fitting’ model gives accurate results for the bulk mobil¬ 
ity of a dumbbell-shaped colloid. Part II of our analysis 
is presented in Ref. [7] and extends these conclusions to 
raspberry particles under confinement. We thus demon¬ 
strate that for a wide range of suitably chosen parameters 
our ‘filling -|- fitting’ formalism leads to a substantially 
improved (and acceptable) numerical tolerance in simu¬ 
lating solid objects with respect to that of the traditional 
raspberry model of Refs. [HI ITS]. 

The remainder of this manuscript is structured as fol¬ 
lows. In Section |ll] we describe our simulation methods 
in detail. Section al A| opens with a description of the LB 
and coupling method. Section fll B| introduces our variant 
of the raspberry model for the spherical and dumbbell¬ 
shaped colloids of interest. Sections |II C| and |II D| detail 
the molecular dynamics and LB simulation parameters, 
respectively. Section HE describes the various hydrody¬ 
namic experiments that we performed to determine the 
properties of the raspberry model. In Section m we 
discuss the dimensionless numbers that characterize the 
physics of our systems. We provide a summary of the 
notations used throughout the text in Section IIG to aid 
the reader when going through the manuscript. In Sec¬ 
tion [ITT] we list our main results. We begin by examining 
the properties of the spherical raspberry in a simple-cubic 
lattice in Section [ill A| We continue with the properties 
of two dumbbell-shaped raspberries with different geo¬ 
metric parametrizations in Section [ill B[ The results are 


discussed and related to previous studies in Section IV 


Finally, we give a summary, conclusions, and an outlook 
in Section El 


II. METHODS 

In this section, we outline the approach used to de¬ 
termine the hydrodynamic properties of a colloid. We 
have split this into subsections detailing the properties 
of the lattice-Boltzmann method, the construction of the 
raspberry model, the molecular dynamics and lattice- 
Boltzmann parameters used, the hydrodynamic exper¬ 
iments performed to extract the mobility of the rasp¬ 
berry, the dimensionless numbers that characterize the 
fluid, and a reference list of the input parameters and 
measured quantities. 


A. The Lattice-Boltzmann Method 

In this section, we briefly outline the major features of 
the lattice-Boltzmann method and viscous particle cou¬ 


pling to put our work into context. We refer the inter¬ 
ested reader to, for instance. Ref. [12] for a more in-depth 
treatment. 

The LB method is a numerical simulation technique 
to solve the Boltzmann transport equation. |12| In its 
simplest form the Boltzmann equation can be written as 

9t/(r,v,t) -b V V/(r,v,t) = C(/(r, v, f)), (1) 

where t denotes time, r the position, and v the velocity; 
dt indicates a partial derivative with respect to time, • in¬ 
dicates the dot product, and V the gradient with respect 
to position; and /(r, v,t) is a phase-space probability 
distribution function and C{f{r, v, t)) is the collision op¬ 
erator acting on the distribution function, which models 
the probability redistribution caused by particle interac¬ 
tions. 

The lattice-Boltzmann equation is the discretized form 
of Eq. Q , where the particle velocities are restricted to 
only a few values. The LB ‘particles’ can thus only move 
in a finite number of directions, which are chosen to be 
commensurate with a space-filling lattice. When this 
lattice has sufficient symmetry to fulfill mass and mo¬ 
mentum conservation, the discrete LB equation can be 
used to determine fluid flow, without directly solving the 
Stokes or Navier-Stokes equations, as has been shown 
via the Chapman-Enskog expansion. |33| The physical 
quantities that are of interest, such as the mass density, 
velocity, and pressure, can be recovered from the modes 
of the discrete probability distribution. 

Current implementations of the LB method trace their 
roots to the lattice gas automata that were developed 
in the late 1980s. [H] |3S] The traditional LB method 
was formulated by making an assumption for the form 
of the collision operator, the right-hand side of Eq. Q, 
the most well-known being the single-relaxation scheme 
introduced by Bhatnagar, Gross, and Krook. [3^ The LB 
method has significant advantages over traditionally used 
fluid solvers, as the algorithm is completely local, which 
allows for straightforward parallelization. [T2| Moreover, 
the streaming operator (left-hand side of Eq. 0 ) and 
the collision process can be fully decoupled, leading to 
an algorithm that is elegant in its simplicity. 

The LB method can be connected to a Molecular Dy¬ 
namics solver, in order to model the behavior of par¬ 
ticles suspended in a viscous fluid. One method to 
achieve particle-fluid coupling was proposed by Ahlrichs 
and Diinweg. m The fluid is coupled to embedded MD 
beads via a friction force that depends on the difference 
in velocity between the bead and the fluid 

Fd = “Co (up — uj(rp)), (2) 

where is the friction force, Co is the bare friction co¬ 
efficient, Up is the particle’s velocity, and is the fluid 
velocity that is evaluated at the particle’s position rp. 
Here, the particle’s coordinates are interpolated onto the 
lattice using a tri-linear scheme. [T2| The opposite force 
has to be applied to the fluid to ensure momentum con¬ 
servation. This algorithm is used to couple the beads of 
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the raspberry model to the LB fluid that will be described 
in the next section. 


B. The Raspberry Model 

In this manuscript we study the so-called ‘raspberry’ 
model for particle-fluid interactions, m as shown in 
Fig.[Tl This model relies on discretizing the surface of 
a larger colloid into coupling points, which experience 
a friction force related to the relative velocity of the 
fluid and the coupling points as described above. [13] In 
Ref. [m, 100 points were used to approximate a sphere. 
To ensure a reasonably homogeneous surface coverage 
these were connected to each other by finite extensible 
nonlinear elastic (FENE) potentials. The forces acting on 
the surface beads were forwarded to a central Lennard- 
Jones (LJ) MD bead, via the LJ interaction. A model 
similar in spirit to the one proposed by Lobaskin and 
Diinweg was developed by Chatterji and Horbach. (TS] 
In their construction the surface beads were fixed with 
rigid bonds to the central bead and no FENE potential 
was employed for the surface-surface coupling. 


1. The Hollow Raspberry 

For the construction of the raspberry model in this pa¬ 
per, we combined the approaches of Refs. [Ml [15]. To ho¬ 
mogeneously arrange the MD beads in a spherical shell 
of radius R, we used a separate MD simulation. We 
placed N > [47rMD beads in a cubic simulation 
box with edge length L, LB lattice spacing a, and peri¬ 
odic boundary conditions. The number of MD beads was 
chosen such that on average there is at least one particle 
per lattice site for the LB simulation. To force the beads 
onto a spherical shell we employed a shifted harmonic 
bond potential around the center of the box, rp, which 
will become the center of the raspberry particle that we 
are creating. This potential has the form 

Vharm(r) = ^A:(|r - rp| - i?)^ (3) 

where r is a point in space and K is the spring con¬ 
stant. To ensure that the beads do not overlap and to 
homogenize the surface density, we endowed them with a 
repulsive Weeks-Chandler-Anderson (WCA) interaction 
potential 



where a is the MD base unit of length and is thus equal 
to the bead diameter. 

The MD beads were thermalized using a Langevin 
thermostat with ‘temperature’ l.Oe and friction coeffi¬ 
cient F = I.OmoT”^. Here, e is the MD base unit of 


energy and corresponds to Ik^T, where fee is the Boltz¬ 
mann constant and T is the temperature, r is the MD 
base unit of time, and ttiq the MD base unit of mass 
(mo = T^etT“^). The MD beads were each given a mass 
Img. By geometrically increasing the spring constant 
from K — l.0ea~^ to K = 3,000ecr“^ the MD beads are 
forced onto the spherical shell described by the potential 
in Eq. (|^. We increased K to its final value oi K = 3,000e 
over 100,000 integration steps of length At = O.OOSr. 
These simulations were performed using the MD software 
package ESPResSo. [UliH] Finally, small deviations of 
the MD beads’ radial position with respect to the desired 
distance R were removed by adjusting their radial posi¬ 
tion. The configuration was then ‘frozen in’ by connect¬ 
ing all beads to a central bead via rigid bonds (virtual 
sites). [IH] 

To test the quality of the result, the raspberry was 
checked for large holes in the surface coverage by applying 
a ‘shotgun’ algorithm. We randomly picked 50,000 points 
on the surface of the sphere and calculate the distances to 
the nearest surface bead. We arrived at the distribution 
of MD beads that we used throughout our simulations, 
by repeating this procedure with different initial config¬ 
urations and particle numbers, until we found a system 
for which the maximum hole size was roughly I.Ocr (bead 
diameter). The outcome for a sphere of radius R = 3.Oct 
is shown in the right-hand side of Fig. [l] Here, 202 sur¬ 
face beads were used to obtain a maximal hole diameter 
of I.Ict. In total five variants of the hollow raspberry 
were considered, with radii R = 2.Oct, 2.5ct, 3.Oct, 4.Oct, 
and 5 .Oct and N = 89, 139, 202, 441, and 593, respec¬ 
tively. We also considered a ‘dense shell raspberry’, with 
R = 3.0(7 and N = 924 beads in the shell, which will 
be discussed further in Section fill A 21 Unless otherwise 
specified, whenever we use the term ‘hollow raspberry’ in 
this document and Part H, [7] we refer to the raspberry 
with N = 202 surface beads and radius R = 3.Oct. 

Finally, we should mention that there is an alternative 
method of positioning the coupling points on the shell. 
For the harmonic potential in Eq. ([^ and WCA inter¬ 
actions between the MD beads, a conjugate-gradient de¬ 
scent method can be used to generate a surface coverage 
with minimal defects. [H] 


2. Filling the Raspberry 


We ‘fill’ the hollow-shell raspberry particle by adding 
coupling points to the interior, as outlined in detail be¬ 
low. We first formed a hollow raspberry according to 
the recipe in Section HBl Next, we added N' > 
|’47r(i? — ct/2)^/3o^] beads to the interior of the shell, 
which interact with each other and the shell MD beads 
via the WCA potential of Eq. Q. The force between the 
internal beads themselves was initially capped to l.Oe/cr 
to prevent numerical instabilities. The system was al¬ 
lowed to evolve by making use of a Langevin thermostat 
(/cbT = 1-Oe, F = I.Otoqt”^). The simulation consisted 
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of over 50,000 time steps of length At = O.OOSr. During 
this time the capping value was slowly raised to lOOe. 
This generally resulted in a random configuration with a 
homogeneous distribution of MD beads within the rasp¬ 
berry. These beads were subsequently frozen in place by 
adding rigid (virtual) bonds to the central MD bead. 

We investigated several values of N' and the role of the 
MD beads’ distribution on the model’s ability to repro¬ 
duce the result of Stokes’ equation, see Section |m A3l 
We settled upon a value of N' = 722, resulting in a 
total of Wot = -I- W -I- 1 = 925 MD beads for the 

so-called ‘filled raspberry’ of radius R — 3.Oct. This re¬ 
sult is shown in the left-hand side of Fig. Note that 
we used exactly the same hollow shell to construct our 
filled variant. For the raspberries with radius R = 2.Oct, 
2.5ct, 4.Oct, and 5 .Oct, we used W = 154, 143, 967, and 
1,323 internal coupling points, respectively. To study the 
improvement of the coupling on the internal filling fac¬ 
tor, we also considered several other values of N' for the 
R = 3.Oct raspberry, namely: W = 50, 100, 200, and 400. 
Unless otherwise specified, we refer to the R — 3.Oct and 
N' = 722 model as the ‘filled raspberry,’ both here and 
in Part II. [7] 

Finally, it should be remarked that in the hydro- 
dynamic simulations utilizing the raspberry model, all 
WCA interactions were switched off and only the rigid 
(virtual) bonds remained. This eliminated any non¬ 
hydrodynamic interactions between the raspberry and its 
images in our simulations with periodic boundary condi¬ 
tions for small box lengths {L « 2R). 


3. Constructing a Dumbbell Raspberry 


A dumbbell-shaped raspberry model (filled or hollow) 
is constructed using a procedure that is analogous to the 
one given in Sections IIB 1 and IIB 2[ Instead of a cen¬ 
tral harmonic potential, we used two harmonic potentials 
centered on rp = (0,0, —d/2) and r'p = (0, 0, d/2), with d 
the distance between the sphere centers of the dumbbell 
(the total length of the dumbbell is d -I- 2R). In addi¬ 
tion, a WCA potential had to be added to prevent beads 
from accumulating in the neck of the dumbbell - the re¬ 
gion where the two dumbbell spheres overlap, if d < 2R. 
To accomplish this, we used a WCA potential between 
the center of the dumbbell, located at (0,0,0), and the 
surface MD beads. This potential had the following form 


Uaeck — 




> 21W 


(5) 


where w is the width of the neck and is given by 




( 6 ) 


After letting the MD beads become trapped in the 
dumbbell shell, in the same manner as for the spherical 


shell, they were connected via rigid bonds to a particle 
at the geometric center of the dumbbell. The dumbbell 
may be filled with N' additional beads using the pro¬ 
cedure outlined in Section [II B2| In this paper, we con¬ 
sider two dumbbell-shaped raspberry particles - one with 
d = 5.0(7 and one with d = 7 .Oct; for both the individual 
sphere radius is i? = 3 .Oct - corresponding to a partially 
overlapping configuration and one with the spheres just 
touching, respectively; see Fig. We used {N = 416, 
N' = 598) for d = 5.0 ct and {N = 502, N' = 404) for 
d = 7 .Oct, respectively, to ensure a homogeneous surface 
distribution and filling of the volume. 



FIG. 2. (color online) Representation of the raspberry dumb¬ 
bells used in our simulations, touching (left) and overlapping 
(right), respectively. The distance between the centers of the 
spheres (each R — 3.0(7 in size) is indicated using the arrows. 
Note that we used the effective MD bead diameter of approx¬ 
imately 1(7 to visualize our result. 


C. Molecular Dynamics Parameters 

Once we had constructed the raspberries, we could use 
them in our LB simulations. The raspberry particles 
were allowed to freely move and rotate, unless other¬ 
wise specihed. All the forces acting on the MD beads 
are transferred to the central bead via the virtual sites 
(rigid bonds). To stabilize the simulation for the bare 
friction coefficients used, we set the (bare) mass and ro¬ 
tational inertia of the raspberry; these quantities should 
not be confused with the virtual mass of the body in a 
fluid, see, e.g., Ref. [SD] for the definition. The mass and 
rotational inertia are based on the particle’s dimensions 
and the fluid mass density, which we denote by p and set 
to p = I.OtoqCT”^. We thus assume that the raspberry 
particle has the same density as the surrounding fluid. 
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For the spheres with radii R = 2.0a, 2.5a, 3.Oct, 4.Oct, used 

and 5 .Oct the mass we used, was m = (4/3)7r/9i?^ « 

33.5mo, 65.5mo, 113mo, 268mo, and 524mo, respec¬ 
tively. The inertia tensor is a diagonal tensor with iden¬ 
tical entries of / = (8/15)7rpii® « 53.6 tooct^, 164moCT^, 

407tooct^, 1.72 • lO^moCT^, and 5.23 • IO^tooct^ for these 
radii, respectively. For the two dumbbell raspberries, we 

The dumbbell’s rotational inertia tensor is diagonal, but 
the entries are not identical. Let /y denote the moment 
for rotation about the main axis of the dumbbell and I± 
the moment for rotation about a central axis perpendic¬ 
ular to the main axis. We may then write 


np ( + dR^ — ) 0 < d < 27? 

12 ; (7) 

8 q 

-irpR^ d > 27? 


'X = 


III = 


7 = 


TTp ( 


np ( 


TTp ( 



I± 

0 0 

0 

/x 0 

0 

0 7|| 


-dR^ + -d^R^ 

4 3 

2 

r 

ldR‘^ - T^-d^7?2 
2 12 


j 3 p2 


1 


24 


960 


i2 n3 


160 


0 < d < 27? 

d> 27? 

0 < d < 27? 
d> 27? 


( 8 ) 


(9) 


( 10 ) 


where the long axis is assumed to be aligned with the 
x-axis. This gives us the following for the d = 5.0 ct 
dumbbell: m « 221mo, I± « 2.23 • lO^moCT^, and 
7|| « 810 tooct^. Whereas for the d = 7 .Oct dumbbell 
we obtain: m « 226too, I± « 3.59 • lO^moCT^, and 
7|| « 814moCT^. 


D. Lattice-Boltzmann Parameters 


The raspberry particles were coupled to an LB fluid 
using the coupling described in Section IIA We did not 
employ the coupling scheme of Refs. mm , since our 
method turned out to work sufficient well for the long¬ 
time properties without modifications to the Ahlrichs and 
Diinweg LB coupling. We used a graphics processing unit 
(GPU) based LB solver, [5T] which is attached to the MD 
software ESPResSo. [13148] The GPU variant of LB im¬ 
plemented in ESPResSo utilizes a D3QI9 lattice and a 
fluctuating multi-relaxation time (MRT) collision opera¬ 
tor. [55] This fluctuating LB model was introduced first 
by Adhikari et al. |53| and later validated by Diinweg et 
al. IMIISSI 


To keep our result as general as possible, we set the 
density of the fluid to p = lmQa~^, the lattice spacing 
to Ict, the time step to At = 0.005t, the (kinematic) 
viscosity to v = la^T~^, the bare particle-fluid friction 
to Co = 25moT“^, and the strength of the fluctuations 


to /cbT = O.Ole, unless otherwise specified. Here, we 
chose neither to optimize our parameters for the most 
accurate reproduction of hydrodynamic interactions, nor 
to match a specific experimental system of interest via 
telescoping. [Ml El] We instead simply chose to use pa¬ 
rameters that are in the regime, where LB reproduces 
hydrodynamic effects for colloids reasonably well and is 
sufficiently stable to use the (float-precision) GPU algo¬ 
rithm, as we will further discuss in Section [lIFj The low 
amplitude of the fluctuations in the thermalized LB is to 
allow averaging over long times without noise dominating 
our results. This will become more clear when we discuss 
these results and prove the importance for the thermal 
averaging performed in Part II [?]■ 


E. Hydrodynamic Experiments 


To assess the quality of the raspberry approximation 
in modeling the hydrodynamic properties of a colloid we 
performed several experiments. We use the term ‘quies¬ 
cent’ to describe an un-thermalized (non-fluctuating, de¬ 
terministic) LB fluid. Below we specify the experiments 
performed for raspberry particles in a simple cubic lat¬ 
tice, i.e., a cubic simulation box of length L with periodic 
boundary conditions. In all experiments the particle was 
initialized in the center of the box. 
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FIG. 3. (color online) Visualization of the various hydrody¬ 
namic experiments carried out in a cubic box of length L 
with periodic boundary conditions. A two-dimensional (2D) 
representation is given here. The blue arrows and symbols 
denote quantities applied to the fluid and raspberry, the red 
arrows and symbols indicate measured quantities. The black 
arrows indicate a thermalized fluid. We refer to the text for 
a description of the experiments, as well as the applied and 
measured quantities. 


• A force experiment in a quiescent fluid, see 
Fig. §a). A constant force F was applied to the 
particle (typically along one of the box axes) and a 


counter force of —F/L^ was applied homogeneously 
to the fluid to ensure that there is no motion of 
the center of mass, i.e., no net transfer of momen¬ 
tum to the system. Not applying this counter force 
would result in an acceleration of the colloid via the 
fluid flow that builds up, as momentum is continu¬ 
ously pumped into the system. The resulting time- 
dependent velocity v(t) and steady-state (terminal) 
velocity Vj were measured and used to determine 
the translational mobility 



( 11 ) 


To establish the steady-state velocity it proved nec¬ 
essary to average over several (very small) oscilla¬ 
tions in the velocity v(t) that are caused by lattice- 
discretization artifacts. 


• A force experiment in a thermalized fluid, see 
Fig. §b). The set-up is the same as for the ex¬ 
periment in Fig. [^a). However, the system was 
first equilibrated until a steady-state emerged and 
the particle fluctuated with the proper thermal ve¬ 
locity distribution. During the production run, v(t) 
was averaged to determine the average steady-state 
velocity Vt = (v(t)), where (•) denotes the time 
average. This allowed us to determine the time- 
averaged translational mobility 



( 12 ) 


• A torque experiment in a quiescent fluid, see 
Fig-iu: c). A constant torque T was applied to the 
particle (typically along one of the box axes). The 
resulting time-dependent angular velocity ojif) and 
steady-state angular velocity were measured and 
used to determine the rotational mobility 


, R _ VM 

Ml i^i 


(13) 


There is no need to apply a ‘back torque density’ to 
the fluid in this experiment, as the periodic bound¬ 
ary conditions do not allow the fluid to develop a 
net rotation. Here, averaging of the oscillations in 
a;(t) due to lattice artifacts also proved necessary. 


• A torque experiment in a thermalized fluid, see 
Fig. id). The set-up is the same as for the exper¬ 
iment in Fig. ic). However, the system was first 
equilibrated until a steady-state emerged and the 
particle fluctuated with the proper thermal distri¬ 
bution. During the production run, <jj{t) was aver¬ 
aged to determine the average steady-state angular 
velocity = (w(t)). This allowed us to determine 
the time-averaged rotational mobility 


-R _ 

Ml irp. 


( 14 ) 
















• A velocity experiment in a quiescent fluid, see 
Fig-iu: e). An instantaneous velocity Vq was im¬ 
parted onto the particle at t = 0 and an instanta¬ 
neous counter velocity of —Vq/L^ was applied ho¬ 
mogeneously (at the same time) to the fluid to en¬ 
sure zero net motion of the system. The resulting 
time-dependent velocity v(t) was measured. This 
quantity can be related to a non-dimensionalized 
velocity auto-correlation function (VACF) 




Vo • v(t) 


|vo| 


(15) 


• An angular velocity experiment in a quiescent fluid, 
see Fig.j^f). An instantaneous angular velocity loq 
was imparted onto the particle at t = 0. The re¬ 
sulting time-dependent angular velocity u}(t) was 
measured. This quantity can be related to a non- 
dimensionalized angular velocity auto-correlation 
function (AVACF) 

asm ^ ( 16 ) 


• An auto-correlation experiment in a thermalized 
fluid, see Fig. |^g). The system was equilibrated 
until the particle fluctuated with the proper ther¬ 
mal distribution. The (A)VACF and the mean 
square displacement (MSD) were measured using 
the multiple-tau correlator in ESPRcsSo. [SB] For 
the (A)VACF the (angular) velocity in the co¬ 
rotating frame was averaged. The C^it) = (v(t) • 
v(t -f r)) and C^it) = (w(t) ■ o}{t -b t)) that fol¬ 
low from the ther mal expe rim ents differ slightly 
from those in Eqs. |l^ and (16), because C'’^(0) = 
SfcBT/m and (7^(0) = as a consequence of 

the equipartition theorem. This allows us to com¬ 
pute the translational and rotational mobility, re¬ 
spectively, via the Green-Kubo relation 


Pl = 


1 




(t)dt, 


(17) 


where the factor 1/3 is used for spherical parti¬ 
cles only and X can be either T or R. [59] The 
relations for anisotropic particles are similar, but 
slightly more involved, since the dot product for 
the (A)VACF is replaced by the dyadic product. 


F. Dimensionless Numbers for the Fluid Properties 


with V the maximum/typical velocity. This implies that 
we can compare it to analytic and numerical results 
obtained by solving the Stokes equations, as will be 
discussed further in Section EH For the colloid radius 
R — 3.Oct and our value of the kinematic viscosity, we 
ensured that the maximum particle velocity remained un¬ 
der O.ISctt”^, for which Re^ < 0.5. However, this value 
was only attained in the velocity and auto-correlation 
experiments for the first time step. For t ^ It and in 
the other experiments, the Reynolds number remained 
smaller than 0.1. Similarly, the rotational Reynolds num¬ 
ber 


Re^ 


ujR^ 

1/ 


(19) 


with uj the maximum angular velocity, remained small: 
Re^ < 0.7, but typically smaller than 0.1. For the 
other radii that we considered, the maximum value of 
the Reynolds number was kept smaller. 

There are a number of relevant parameters to describe 
the hydrodynamic properties of our system. For the ther¬ 
malized LB fluid, we can define the Peclet and Schmidt 
number of the particle, and the Boltzmann number of the 
fluid. In both quiescent and thermalized LB fluids, we 
can determine the Mach number. Finally, the coupling 
of the raspberry particles to the fluid can be described 
by the Immersion number and the Screening ratio. We 
will determine these numbers next. 

The translational and rotational Peclet numbers are 
defined as 


Pe^ 

Pe^ 


vR 

^0 


nfl ■ 
^0 


( 20 ) 

( 21 ) 


For the thermalized force and torque experiments, see 
Figs. [^b,d), we obtain values of Pe^ « 1.5 • 10^ and 
Pe^ « 5.0 • 10^ for the R = 3.Oct raspberry. We did not 
carry out similar thermalized experiments for R ^ 3.Oct. 
If we use the thermal velocity for the auto-correlation 
experiment, see Fig. [^g), then we arrive at Pe^ « 2.5 • 
10^ and Pe^ ~ 5.0 • 10^. The large value of the Peclet 
number indicates that our results are in a regime that 
is dominated by advective flow, rather than by diffusion. 
That is, our thermalized results can be readily compared 
to those of the quiescent (deterministic) experiments. 

The Schmidt number of the particles measures the rel¬ 
ative importance of diffusive and advective transport and 
is defined as 


Sc = 



( 22 ) 


In the above experiments, care was taken to en¬ 
sure that the particle remained in the low translational 
Reynolds number regime 

Re^ = — <C 1, (18) 

V 


where v is the kinematic viscosity, as before. We obtain 
Sc = 5.6 • 10^ for the R = 3.Oct colloid and a thermal 
fluctuation strength of k-^T = O.Ole. This value of the 
Schmidt number is quite high, compared to the typical 
value of Sc « 10 in LB simulations. However, it is nec¬ 
essary to use such high numbers to access the regime in 
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which the momentum diffusion in the fluid dominates the 
diffusive transport of the particles. This allows for the 
accurate measurement of hydrodynamic interactions in 
confining geometries, see Ref. [7]. 

The Boltzmann number Bo of the LB fluid, which in¬ 
dicates the level of coarse graining, is defined as 


Bo = - 


{{vf) - {Vj? 

{v^) 


> 1/2 


(23) 


where Vi is the speed of the LB fluid at a given node in the 
fluid. [12] By averaging over 10^ LB nodes for the param¬ 
eters that we used, we obtain Bo = 0.82. For Bo = 1 (the 
maximum value) the model is fully microscopic, whereas 
for Bo = 0 the model is entirely deterministic. For this 
value of the Boltzmann number we are in an intermediate 
regime, with a limited level of coarse graining. 

The Mach number of the LB fluid is the ratio of the 
particle velocity to the speed of sound and is given by 

Ma=—, (24) 

Vs 

where Vs is the speed of sound in LB 


R 

Xtot 

Co 

Kx 

Rl 

■nR 

XXp 

Rl/R^ 

Filled 

5.0(7 

1917 

25moT~^ 

32.3 

4.84(7 

4.85(7 

0.9986 

4.0(7 

1409 

25moT~^ 

30.9 

3.86(7 

3.87(7 

0.9985 

3.0(7 

925 

35moT~^ 

31.0 

2.90(7 

2.90(7 

0.9985 

3.0(7 

925 

25moT~^ 

28.9 

2.89(7 

2.90(7 

0.9983 

3.0(7 

925 

10moT~^ 

22.3 

2.86(7 

2.87(7 

0.9971 

3.0(7 

603 

25moT~^ 

23.4 

2.86(7 

2.87(7 

0.9973 

3.0(7 

403 

25moT~^ 

19.1 

2.83(7 

2.84(7 

0.9960 

3.0(7 

303 

25moT~^ 

16.6 

2.80(7 

2.82(7 

0.9948 

3.0(7 

253 

25moT~^ 

15.1 

2.78(7 

2.80(7 

0.9938 

2.5(7 

248 

25moT~^ 

16.4 

2.34(7 

2.35(7 

0.9947 

2.0(7 

245 

25moT~^ 

18.2 

1.88(7 

1.89(7 

0.9957 


Hollow 


5.0(7 

594 

25moT~^ 

18.0 

4.93(7 

4.95(7 

0.9953 

4.0(7 

442 

25moT~^ 

17.3 

3.94(7 

3.96(7 

0.9950 

3.0(7 

203 

35moT~^ 

14.5 

2.94(7 

2.96(7 

0.9928 

3.0(7 

925 

25moT~^ 

28.9 

2.98(7 

2.99(7 

0.9982 

3.0(7 

203 

25moT~^ 

13.5 

2.93(7 

2.95(7 

0.9918 

3.0(7 

203 

10moT~^ 

10.4 

2.88(7 

2.92(7 

0.9861 

2.5(7 

140 

25moT~^ 

12.3 

2.42(7 

2.45(7 

0.9900 

2.0(7 

90 

25moT~^ 

11.0 

1.93(7 

1.95(7 

0.9876 


Vs = 



(25) 


with a the lattice spacing and At the time step. The 
prefactor depends on the shape and dimensionality of 
the grid (the prefactor for a D2Q9 grid is the same inci¬ 
dentally). For our parameters we obtain Ma~ 9.0 • 10“^ 
for the thermalized force experiment of Fig. ib) and 
Ma ~ 1.5 • 10“"^ if we take v to be the thermal velocity 
of a i? = 3 .Oct colloid at kuT = O.Ole. For all radii we 
obtain Ma < 10“^. 

The immersion number, which measures the relative 
density of the MD beads, is defined as 


In = 



(26) 


where TOq is the mass of a single MD bead. It should 
be noted that the individual (virtual) MD beads, which 
make up the raspberry, have unit mass (mbead = l.Omo). 
Only the central bead, to which the other beads couple 
and which holds the properties of the entire colloid, has 
a different mass on the MD level. For our choice of pa¬ 
rameters we obtain In = 0.5, which corresponds to a 
neutrally buoyant object. 

The screening ratio |5UMT] for a filled sphere is a mea¬ 
sure for the porosity of an object and it is given by 


TABLE I. The screening-ratio related properties for the vari¬ 
ous raspberry particles studied in this manuscript. From left 
to right, the columns give: the imposed radius R, the total 
number of coupling points in the raspberry Ntot, the value of 
the bare friction coefficient (^o, the screening ratio fcx (‘X’ is 
either ‘F’ or ‘H’), the translational hydrodynamic radius for 
the porous sphere Rp, |391140| the rotational hydrodynamic 
radius for the porous sphere Rp, El and finally the ratio of 
these radii. 


is the density of coupling beads in the sphere (assuming 
a uniform distribution - which is an acceptable approxi¬ 
mation for our fillings), (^eff is the effective single particle 
fluid-coupling [12] 


Ceff 


(m-V 

VCo w ) 


(29) 


with ^0 the imposed LB fluid friction, 0.04 a prefactor for 
the D3Q19 grid, the value of which we measured, and a 
is the LB grid spacing. For the various particles that we 
used the screening ratio is given in Table [l| For a hollow 
sphere [39| |40| the screening ratio is defined by 


kh = 



(30) 


II 

05 

(27) "'‘'™ 


V p 

Ntot 

where ry is the dynamic viscosity, and 




(31) 


Pf 


3Atot 


is the density of beads on the surface of the sphere as- 

(28) 

suming a uniform distribution. 
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The screening radio gives insight into the match be¬ 
tween the translational and rotational hydrodynamic ra¬ 
dius of the particle, as predicted by porous particle the¬ 
ory. For the filled sphere the following hydrody¬ 

namic radii are expected 


RI = R- 


^ tanh(KF)^ 




1 + 


2k| 


1 - 


tanh(KF^ 


Kf 


= i? ( 1 - 3 


COth(KF) 


Kp 


1/3 


(32) 


(33) 


where Rp |39j is the translational hydrodynamic radius 
and Rp [JT] is the rotational one. For the hollow shell 
the expressions are [MMT] 


Rp = i?-g- 


2 +«H 


R^ = R 


1/3 


(34) 

(35) 


The values of these radii are given in Table |T] It should 
be noted that the ratio of the two radii (translational 
and rotational) is only equal when kx t oo (‘X’ is either 
‘F’ or ‘H’). For kx i 0 the translational hydrodynamic 
radius is much smaller than the rotational one. We have 
added the values of these ratios to Table|T| In our analysis 
and discussion, see Sections |III A 3| and |W] we relate the 
insights of porous-particle theory to our results for the 
raspberry particle. 


the angular velocity auto-correlation func¬ 
tion (AVACF) for rotation in a cubic box of 
length L with periodic boundary conditions, see 
Figs. [^f,g) and Eq. ([^). 


• the time-dependent translational mobility in 

a cubic box of length L with periodic boundary con¬ 
ditions, see Fig. HJa). When the time dependence 
is dropped, the limit t f oo has been taken. 


• /(^(t), the time-dependent rotational mobility in a 
cubic box of length L with periodic boundary con¬ 
ditions, see Fig. [^c). When the time dependence 
is dropped, the limit t f oo has been taken. 


• the time-averaged translational mobility re¬ 
sulting from the thermal force experiment, see 
Fig-db). 

• pL^, the time-averaged rotational mobility resulting 
from the thermal torque experiment, see Fig. j^d). 

• /(q , the bulk translational mobility, which follows 
from the limit L f oo of 

• the bulk rotational mobility, which follows from 
the limit L f oo of 


• /, the fractional deviation between two results. 


III. RESULTS 

In this section, we discuss the results that we obtained 
by performing the simulations and numerical calculations 
outlined in Section [H] We have split this section into two 
parts: one for the sphere and one for the dumbbell. These 
parts are further subdivided according to the nature of 
the experiments. 


G. Notations Used throughout this Manuscript 


In this section, we summarize the notations used in this 
manuscript. This will aid the reader in going through the 
text, as many of the notations are necessarily similar. 


• L, the box length of the cubic box with periodic 
boundary conditions. 

• the effective hydrodynamic radius obtained by 
extrapolating translational mobility measurements, 
see Figs. da,b,e,g), for the limit of box length 
L t oo. The subscript h is used to differentiate Rh 
from the bead-to-center distance of the raspberry’s 
coupling points R. 


• RJ^, the effective hydrodynamic radius obtained by 
extrapolating rotational mobility measurements, 
see Figs.[^c,f,g), for the limit of box length L f oo. 


Clit), the velocity auto-correlation function 
(VACF) for translational movement in a cubic box 
of length L with periodic boundary conditions, see 
Figs. |3];e,g) and Eq. ([^|. 


A. Sphere iu a Simple Cubic Crystal 

1. The (Angular) Velocity Auto-Correlation Function 


Using the (quiescent) velocity and (thermalized) auto¬ 
correlation experiments discussed in Section |IIE[ see 
Eigs. [^e,g), we established the VACE for a filled rasp¬ 
berry sphere in a cubic box of length L — 100.Oct. The 
results are shown in Fig. From Figs. Qa,b,c) we ob¬ 
serve several decay regimes that are typical for the LB 
simulations of the raspberry particle. In the following 
they will be described in more detail. 


Decay Regimes 


(I) At short times there is an unphysical-coupling 
regime, see Fig. Sa), in which the VACF decays expo¬ 
nentially according to 


cm 

cm) 


= exp 




(36) 
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FIG. 4. (color online) The velocity auto-correlation fnnction (VACF) Cj (t) as a fnnction of the time t expressed in the MD time 
unit T. The graphs show results for a filled raspberry of radius R = S.Ocr in a box of length L — 100.Oct, with LB parameters as 
given in the text, (a) The initial decay of the VACF. The red sqnares with error bars show the result for a thermalized LB, the 
blue solid curve gives the result of a quiescent experiment, and the green dashed line shows the predicted unphysical-coupling 
decay inherent to LB. (b) Log-linear plot of the initial and intermediate decay of the VACF. The dashed orange curve gives the 
expected Stokes’ decay. An arrow indicates the position where the correspondence is reasonable. The vertical lines indicate the 
time for sound waves to propagate through the system over certain lengths: one lattice spacing {tc, dashed black), roughly the 
inter-bead separation; the raspberry’s hydrodynamic diameter (t 2 H, dashed gray); and the box length {II, dotted gray), (c) 
Log-log plot of the long-time decay. The magenta line shows the power-law decay. The unphysical coupling, fitted Stokes’, and 
final exponential-decay (purple dashed, indicated by an arrow) curves are shown for completeness. The two brown vertical lines 
indicate the time needed for viscous dissipation over certain lengths: the radius of the sphere (t/j 2 , dashed), and the product of 
the sphere radius R with the box length L {Irl, short dashes), (d) The time-dependent Green-Kubo value of the translational 
mobility obtained from the quiescent (blue solid curve) and thermalized (red squares with error bars) LB result. The 

solid cyan line shows the result of a quiescent force experiment {fjJ,, derived from the terminal velocity vt), while the dashed 
cyan line shows the result of a thermal-averaged force experiment {Rr, from the time-averaged terminal velocity vt)- 


with Ntot the total number of beads, Co the bare friction 
coefficient, and m the particle’s (bare) mass. [2] The ex¬ 
istence of this regime can be attributed to the fluid not 
co-moving with the velocity of the beads (raspberry par¬ 
ticle). That is, the MD beads interact with the stationary 


fluid only through a regular Langevin-type friction - the 
velocity of the fluid is essentially zero during these time 
steps. 


The expected (unphysical) decay of Eq. (36) is indi¬ 
cated in Fig. [^a) and matches reasonably well with the 
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observed initial decay. However, the result deviates even 
in the first and second time step, signifying the onset 
of proper coupling. This is in agreement with the re¬ 
cent observations in the MPCD simulations of Ref. m, 
where this deviation from the expected unphysical decay 
was also attributed to the onset of hydrodynamic corre¬ 
lations. Finally, note that there is a small deviation be¬ 
tween the thermalized LB result and the quiescent VACF 
when t > O.OSr, to which we will come back later. 

From the above it is thus clear that the no-slip bound¬ 
ary condition at the surface of the raspberry is violated 
at short times, even taking the hnite compressibility of 
the LB fluid into account. Moreover, the expected decay 
for a porous colloid [60] is not captured by the raspberry 
with the Ahlrichs and Diinweg coupling. [T3] This is a 
problem inherent to the LB method. HKis] The mod- 
ihed coupling scheme by Mackay et al. |36| purportedly 
remedies this problem, we will come back to this in Sec¬ 
tion |TV| 

(11) The Decay. At intermediate times there is a 
regime, in which the VACF decays exponentially accord¬ 
ing to Stokes’ prediction 


cl{t) 

cm ^ 



dirriR \ 
m / 


(37) 


Here, we used the proportionality symbol, since the un¬ 
physical initial decay makes it impossible to establish an 
analytic prefactor for the onset of this regime in fluid- 
particle coupling. The regime appears because the hydro- 
dynamic coupling between the raspberry particle and the 
surrounding fluid is now fully established. m It should 
be noted that in Ref. mi the mass in the denominator 
was taken to be m*, where m* is the ‘virtual’ mass. [50] 
This virtual mass is the particle mass m plus half of the 
displaced fluid mass; m* = in our case. We will 

come back to this shortly. 

The applicability of Stokes’ prediction for our numeri¬ 
cal results can be seen in Fig. |^b), where a Stokes-type 
decay has been htted to our data. The agreement is not 
very convincing. The curve does not match the Stokes’ 
trend well. However, the agreement between the bare- 
mass prediction of Eq. (371 is superior to the one in 
which the virtual mass is used (not shown here). The 
latter type of decay was originally suggested by Lobaskin 
and Diinweg. m The superiority of the bare-mass re¬ 
sult could be reasonable since Felderhof [SD] has shown 
that for a porous sphere the TO*-related decay regime is 
absent in the high-frequency limit. Unfortunately, it is 
unclear whether our simulations are sufficiently close to 
this limit. In addition, in the limit where the viscous 
coupling constant goes to infinity before the frequency, 
the virtual mass decay is present. [60] The fact that the 
high-frequency porous sphere solution of Felderhof [50] 
does not match better in the Stokes-type regime, makes 
for a slightly academic discussion, since such comparison 
is hindered by the presence of the unphysical decay. 

Characteristic Times. We have indicated three charac¬ 
teristic times related to sound propagation in the LB in 


Fig. I^b). The three times are = cr/vs, t 2 R = 2R/vs, 
and tR = L/vs, i.e., the time required for sound waves to 
propagate one lattice spacing, the diameter of the rasp¬ 
berry, and the length of the box, respectively. Here, Vs is 
the speed of sound, as defined in Eq. (25). We will now 
discuss the relevance of these times. 

For the filled sphere, in which the MD beads are 
roughly l.Ocr apart, we find possible signatures of the 
propagation of sound between the MD beads, as can be 
inferred from the short-time oscillations. The first dip 
in the VACF roughly coincides with ~ 8.7 • 10 “^t, as 
indicated by the black dashed line in Fig. Qb). These 
oscillations may also be related to the magnitude of the 
effective friction that the added coupling points in the 
interior bring about. At the time it takes sound to prop¬ 
agate the diameter of the sphere {t 2 R « 6.1 • 10“^r), we 
find a small dip in the VACF, see the dashed gray line 
in Fig. fflb). This dip is similar to the one observed in 
Ref. IdoT and is caused by the compressibility of the LB 

fluid, m 

Note that the Stokesian regime of decay appears to be 
delimited by the time it takes sound to travel the distance 
of the box (fi « 0.87r, dotted gray line in Fig. Qb)). 
However, for our specific choice of parameters, this time 
is close to the viscous time it takes momentum to dif¬ 
fuse by one colloidal radius tR 2 = pR'^/r] = 9.Or. This 
viscous time is the relevant time scale for the develop¬ 
ment of hydrodynamic memory effects. [SOI Eg We have 
a stricter separation of sonic and viscous time scales than 
in Refs. HgED], i.e., tR 2 /t 2 R » 1. Therefore, our results 
do not display sound undulations (back tracking) in the 
long-time power-law regime. 

(Ill) After a sufficiently long time, the hydrodynamic 
interactions with the surrounding fluid result in a per¬ 
sistence of the velocity (non-exponential decay) as the 
vorticity diffuses away from the particle. These hydro- 
dynamic memory effects lead to an algebraic decay of the 
(A)VACF; the so-called ‘long-time tail’. [5T] This decay 
has the following form 


= ^"iVp(7r77t) 3 / 277 ( 7 ?, L); (38) 

= 7Tl^{Anr,t)-m{R, L), (39) 

for the translational and rotational motion, respec¬ 
tively. [50l [ 59 I ESI [63] N.B. These are the 3D auto¬ 
correlation functions, which are normalized. This 
was unclear in our J. Chem. Phys. publication. 

Here, 77(7?, L) is the Hasimoto scaling expression [5] 


77(7?,L) = 1 - 2.837 (^1^(^1^ . (40) 

Figure |^c) shows the power-law decay for the transla¬ 
tional motion more clearly. The correspondence with the 
quiescent data is excellent, we obtain a match for both 
the prefactor and exponent via a fitting procedure that is 
within 1% of the theoretical prediction. Note that within 
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the error bar, which gives the standard error, the decay 
is captured by the thermalized result. The thermal data 
shows correspondence within the error bar, however the 
error bars are substantial in this regime; it was the best 
that could be achieved within a reasonable time frame 
for our choice of parameters. Only for L 3> 30i? is the 
power-law decay more pronounced. However, larger box 
sizes require even longer sampling. Our result is similar 
to that of Refs. [TillT^ . 

In Fig. I^c) we have indicated two times related to vis¬ 
cous dissipation over certain length-scale combinations: 
tjii = 9.0t (as before) and tuL = pRL/r] = 300.Or. 
These two times roughly indicate the start and end of the 
power-law decay. In Ref. [64] it is suggested that the ex¬ 
ponential decay that follows the power-law decay, will set 
in at 1^2 = /oT ^/77 = 1.0 • lO'^r. However, from Fig. |^c) 
it is clear that this third exponential decay, which will be 
discussed next, sets in far sooner than this. 

(IV) For the quiescent data, there is a third expo¬ 
nential decay in the data when t > Irl, see the purple 
dashed line in Fig. |^c). Analysis shows that this decay 
has a small exponent that depends on the size of the sim¬ 
ulation box. In Ref. |Mj the following form for the decay 
is suggested 


lines in Fig. id), respectively. We arrived at a value 
of = 1.38 • 10“^cr^e“^r“^ for the quiescent data and 
/iD = 1.32 • 10“^CT^e“^T“^ for the thermalized data. The 
results from the VACF and the force experiments cor¬ 
respond within the error, but there is a discrepancy be¬ 
tween the thermal and quiescent data. This deviation 
can be explained by the way these experiments are car¬ 
ried out. The counter-force applied to the fluid also acts 
inside of the particle, effectively modifying the total ap¬ 
plied force, as will be further discussed in Section [TH A 2| 
For completeness we performed quiescent and thermal 
torque experiments, see Figs.|^c,d). From these, we ob¬ 
tained /rD = 9.24- for the quiescent data and 

= 9.25 • for the thermalized data. These 

correspond well within the numerical error. This further 
proves the idea that the mismatch between the and 
/iD can be attributed to an artifact of the experiment. Fi¬ 
nally, it should be noted that a similar mismatch between 
/iD and can occur for small values of the particle’s 
Schmidt number {Sc, Eq. (22)), |HS] which are typical 
for LB. However, as we have shown in Section |HF| our 
particle Scmidt number is quite large {Sc > 10^), which 
should put our system in the regime where the thermal 
and quiescent (deterministic) results correspond well. [65] 


Cj {t) ( 47r^?7 

^o^exp(^-—t 


(41) 


which according to Ref. [M] should set for t > tiji . The 
exponent comes from the smallest positive value over all 
potential wave numbers that fit the geometry of the box. 
We observe that the decay sets in more quickly in our 
simulation, namely around t « Fitting for the value 
of the exponent we obtain 4.7T0“^r“^, whereas the form 


in Eq. (41) yields 3.9-10 There is a difference 


between these two factors of about 20%, but within the 
error the decay is well captured by Eq. (41) - only our 
fit is shown in Fig. [^c). The cause of the early onset of 
the third exponential decay is unclear at this time. 


Thermal versus Quiescent 


Finally, we considered the Green-Kubo relation for the 
VACF by integrating 


J {t')dt', (42) 

for the thermalized data. The expression for the quies¬ 
cent data is similar. Figure|^d) shows the resulting time- 
dependent translational mobility {t). We obtained the 
value of = pL^{t t 00 ) = 1.37 • 10“^(T^e“^T“^ from the 
quiescent data for the box of length L = 100.Ocr. The 
data for the thermalized LB has a slightly lower value 
than the quiescent result, which can in part be attributed 
to the deviation that was already present at short times. 

In addition to determining from the VACFs we per¬ 
formed a quiescent and thermalized force experiment. 
The result is shown using the solid and dashed cyan 


The Hollow Raspberry 


In order to examine the difference between the hollow 
and filled raspberry model, we carried out similar exper¬ 
iments for a hollow-raspberry sphere in a box of length 
L = 80.Oct. We find similar regimes as in Fig. For 
the hollow raspberry there is weaker coupling with the 
fluid. This is caused by the spatial distribution of cou¬ 
pling points and reduced the number of points. This re¬ 
sults in weaker decay of the unphysical-coupling regime, 
which therefore matches the exponential form of Eq. (36) 
more closely. 

Note that the existence of the power-law behavior 
is more convincingly shown by our AVACF data, see 
Fig.gb), as the fitted function and measured decay cor¬ 
respond well over a decade in time. It is unclear whether 
the modified coupling scheme by Mackay et al. [521 shows 
a similar decay. Finally, it should be noted that for the 
hollow raspberry VACF there is the same deviation be¬ 
tween the quiescent and thermalized LB results as shown 
in Fig.|^a). However, the thermal and quiescent data for 
the AVACF match well throughout. 


2. The Influence of Crystal Lattice Spacing 

Thus far, we have examined only the results of the (an¬ 
gular) velocity experiments, and shown that these cor¬ 
respond - at least amongst themselves - to the results 
of force experiments for the same system. Let us now 
consider the effect of the lattice spacing of the simple 
cubic crystal on the hydrodynamic coupling between the 
spheres. This simple cubic geometry is unlike a physical 
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FIG. 5. (color online) Auto-correlation functions for a hollow 
raspberry of radius R = 3.0cr in a box of length L = 80.Oct as 
a function of the time t expressed in the MD time unit r. The 
LB parameters are as given in the text. The use of symbols is 
the same as in Fig. (a) Velocity auto-correlation function 
(VACF) for the thermal (red squares with error bars) 

and quiescent (blue curve) calculations. The unphysical cou¬ 
pling (dashed green), the fitted Stokes (orange dashed), fitted 
final exponential (purple dashed) curves, and power-law decay 
(magenta) are also shown. The two gray dashed vertical lines 
show the time it takes for sound to travel the particle’s hydro- 
dynamic diameter (t 2 H, dashed) and the box length {II, dot¬ 
ted), respectively. The two brown dashed vertical lines show 
the time associated with viscous dissipation (t^ 2 , dashed) and 
{tRL, short dashes), (b) The angular-velocity auto-correlation 
function (AVACF) C^{t) for the same parameters. 


crystal, in the sense that all particles translate and rotate 
in unison; an effect of there being only a single particle 
in a box with periodic boundary conditions. However, 
there is experimental evidence that such systems may be 


realized. |66H68) The uniformity of the periodic structure 
makes the solutions to Stokes’ equations for this geom¬ 
etry analytically tractable. Such calculations were per¬ 
formed, for example, in the work of Hasimoto [5] and of 
Hofman et al. [1] 


Translation of the Crystal 


Figure [^a) shows the change in velocity v{t) during 
a quiescent force experiment (see Fig. |^a)) for a num¬ 
ber of box sizes L using the filled raspberry model. Note 
that for larger L the friction experienced by the parti¬ 
cle is smaller, as the hydrodynamic-interaction with its 
periodic images is reduced. However, the time it takes 
for the stationary state to set in is increased, as it takes 
longer to transfer momentum between the particle and 
its images. From the terminal velocities in the station¬ 
ary state we determined the mobility, by averaging over 
several oscillations due to lattice artifacts. 

In order to establish the mobility at infinite dilution 
(one particle in bulk), we fitted our data using a poly¬ 
nomial of the Hasimoto form: 0 A + BfL -I- C/L^, see 
Eq. (40), in the range where this form is expected to be 
valid (L > 3i?, as we will see later in this section) and 
extrapolated to L f oo. The resulting value for A is the 
bulk translational mobility 


T _ 

Mo = 


dnrjR^ 


T ’ 


(43) 


with the translational hydrodynamic radius. We were 
thus able to determine the extrapolated value and si¬ 
multaneously the effective hydrodynamic radius of our 
raspberry colloid, using Eq. (43). This extrapolation 


refers to the ‘fitting’ part of our ‘filling -|- fitting’ for¬ 
malism. These two parameters and allowed us 
to non-dimensionalize the box length and the measured 
translational mobility, as shown in Fig. i3. 

In Figs. [^b,c) we compare the quality of our result 
for the box-size dependence with the analytic result by 
Hasimoto [8] given in Eq. (40) (dashed red curve) and 
the numerical calculations by Hofman et al. [1] (dashed 
green curve). Figure |^c) shows the fractional deviation 
/ between our data and the two literature results, as 
well as the difference between the Hasimoto (Ha) and 
Hofman et al. (Ho) data. For the data points provided 
by Hofman et al. [1] we used a polynomial fit to represent 
these as a curve. The fit has the following shape 


{R, L) = 1 - 2.807 - -t- 3.437 - 


R 


R 


(44) 


Note that the analytic and numerical expressions of 
Refs. [H [S] correspond well for box sizes greater than 
L « 5.0i?^. That is, within the error expected for the 
fitting procedure that we applied to the data by Hof¬ 
man et al., there is good agreement between their results 
and Hasimoto’s data over this range. The discrepancy 
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FIG. 6. (color online) (a) The velocity v{t) as a function of the time t expressed in the MD time unit r obtained from force 
experiments for a selection of box sizes L. (b) The dependence of the particle’s translational mobility (expressed in terms of 
the bulk translational mobility hq) on the inverse box length 1/L times the hydrodynamic diameter (twice the hydrodynamic 
radius The blue circles show results obtained from the velocity experiment. The red dashed curve shows the analytic 

expression by Hasimoto [S] (Eq. (40l) and the green dashed curve shows a polynomial fit to the numerical data of Hofman et 
at, |4] see Eq. (441. The black crosses give the results of the force experiment, corrected for the counter force applied inside 
of the raspberry, as explained in detail in the main text. The vertical dashed black line indicates the value of L for which the 
spheres are separated by one lattice spacing (L = 2RJ^ + a), (c) Fractional deviation / as a function of jL. The solid blue 

curve shows the difference between the theoretical expressions of Hasimoto and Hofman et at. The red circles and green squares 
indicate the difference between our raspberry experiment and the Hasimoto and Hofman et al. expressions, respectively, (d) 
Fractional deviation / for the force-corrected data. The blue line is the same as in (c). The magenta pluses and black crosses 
give the difference between the force-corrected data and the Hasimoto and Hofman et al. expressions, respectively. The gray 
horizontal line in (c,d) indicates a fractional deviation of 2.5%. 


for smaller box sizes can be explained by the truncation 
of the series expansion in Hasimoto’s work. 

Our raspberry results (Ra) agree reasonably well with 
the data of Hofman et al. over the range L > 5.0i?^, 
but there is also a clear signature of systematic deviation 
present in /. This implies that our data differs substan¬ 
tially from the values of Ref. [1] in the 1/L^ term. A 


similar range of agreement and small-box-size deviation 
can be observed between our data and that of Hasimoto. 
However, in spite of this, our data is much closer to the 
results of Hofman et al. than those of Hasimoto; by al¬ 
most an order of magnitude in / for L | 2i?^. We will 
discuss the origin of this systematic deviation between 
our data and that of Ref. [4] next. 
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Origin of the Discrepancy 

The discrepancy between our data and the result by 
Hofman et al. brings us back to the difference that 
we observed between the VACFs obtained from the ve¬ 
locity and temperature experiments carried out in Sec- 
Fig. ga). Remember that in the quies- 


tion 


IIIAl 


see 


cent experiments a homogeneous and instantaneous ve¬ 
locity has to be applied to the fluid in order to ensure 
zero net movement of the system, see Fig.ge). Similarly, 
for the quiescent force experiment, a constant homoge¬ 
neous force density is applied to the fluid, see Fig. ga). 
Consequently, this velocity and force are also applied di¬ 
rectly to the fluid nodes that are coupled to the rasp¬ 
berry MD beads. The effective force applied to the colloid 
can therefore be calculated by subtracting the integrated 
fluid force-field over the volume of the raspberry. This 
calculation yields 


feff = f 1 - 


47ri?3\ 

) ’ 


(45) 


where f is the force directly applied to the central bead of 
our raspberry construct. Analogously, the counter veloc¬ 
ity affects the time evolution of the VACF. For the ther- 
malized experiments this was not an issue, since counter 
velocities and forces do not need to be applied. These 
counter velocities and forces are therefore a likely can¬ 
didate for the observed discrepancies. This implies that 
the force/velocity experiments are unsuited to analyze 
the hydrodynamic properties of hnite systems in their 
present form. The fact that there is a mismatch between 
the thermal and quiescent results in Fig. ga) is thus not 
an expression of a violation of the equipartition theorem 
or fluctuation dissipation. Nor is it correct to argue that 
this is a consequence of the porosity of the particle. The 
counter-force is only used to counter momentum transfer 
to the periodic system by the force applied to the parti¬ 
cle. The behavior in the limit of the infinite system is, 
however, accurately captured, as the back velocity and 
force vanish. 


We took the effective force of Eq. (451 to determine the 
‘corrected’ value of (Co) using Eq. (11) as a function 


of the box size, see Fig. gb). Note that the correspon¬ 
dence between the result by Hofman et al. and our data is 
thus greatly improved and that the systematic deviation 
is removed for large box sizes, see Fig. gd). Moreover, 
for small box sizes the deviation between our corrected 
result and the literature values is substantially reduced, 
although a systematic difference remains. Within the 
error, the data corresponds much closer to the data by 
Hofman et al. than it does to the Hasimoto result. 

From our corrected data, we estimated the range over 
which the raspberry is able to accurately reproduce hy¬ 
drodynamics interactions (/ < 2.5%) in our system. 
For this particular model we found the criterion to be 
L > 2.8RJ^, which can be extrapolated to other spatial 
arrangements of the colloids. It is likely that this crite¬ 


rion can be extended to smaller boxes, as we will see in 
the following and in Part H. [7]. The normalized results 
for a hollow raspberry lie on top of the filled ones shown 
in Fig. gb) within the error bar (not shown here). How¬ 
ever, the values for the effective hydrodynamic radii 
differ: 3.53 (t and 3.47 ct for the filled and hollow model, 
respectively. 

(a) 


DC o 

DCtl] 

=t 


1 

1 1 

1 1 

1 

1 

1 

0.9 


V 1 

0.8 


R 1 
! 

\ ! 

0.7 


\! - 



rt! 

0.6 


1 

i\ “ 

0.5 

Filled o 

Hollow □ 

Hofman .. 

_1_1_ 

1 \ 

1 \ 

_ ^ 


0.2 


0.4 


0.6 


0.8 


2Rh/L 


(b) 


10"^ 


10" r 


10"" r 


10“^ b- 


10 


-4 


Filled 

Hollow 


o 

□ 


-S 




□ 

o 


o 

oO 

o 

□ o 


_L 


_L 


0.2 


0.4 


0.6 


0.8 


2R^/L 


FIG. 7. (color online) (a) The dependence of the particle’s 
rotational mobility gf; (expressed in terms of the bulk rota¬ 
tional mobility ^q) on the inverse box length times the hydro- 
dynamic diameter (twice the hydrodynamic radius i?^). The 
blue circles show results obtained for the filled raspberry and 
the green squares for the hollow raspberry. The red dashed 
curve shows the expression given by Hofman et al, [4] see 
Eq. (461. The vertical dashed black line indicates the value of 


L, for which the spheres are separated by one lattice spacing, 
(b) Fractional deviation / as a function of 2R^/L. The blue 
circles and green squares indicate the difference between the 
filled and hollow raspberry and analytic expression, respec¬ 
tively. The gray horizontal line indicates a fractional devia¬ 
tion of 2.5%. 
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Rotation in the Crystal 

We continued our verification of the quality of the filled 
and hollow raspberry model by examining hydrodynamic 
coupling between spheres rotating in unison in a cubic 
lattice, as before, see Fig.|^c). Figurej^shows a compar¬ 
ison of our results to the expression given by Hofman et 
al. [1] for the box-size dependence of the rotational mo¬ 
bility The expression provided in Ref. [J] reads 

= 1-4.189 


/ rj>\ 10 

-h 231.858 f - j . (46) 



The procedure used to generate this data is analogous to 
that outlined for the translational experiments. Using 


= 


Sirri (i? 


(47) 


we determined the effective hydrodynamic radius 
from our data. Note that while there is still a systematic 
component to /, see Fig. 00, the agreement between 
our result and literature is excellent for both models. 

This further demonstrates the plausibility of our as¬ 
sertion that the high level of deviation for the transla¬ 
tional mobility is caused by the back-force/velocity that 
is applied homogeneously to the fluid, since a similar cor¬ 
rection is not required for the rotational experiments. 
However, there is a fundamental difference between the 
experiments. The rotational motion exposes the fluid 
to constantly varying coupling points (the MD beads), 
whereas for translational motion the fluid could more 
easily find a pathway of least resistance. This could be 
another source of discrepancies in the translational ex¬ 
periments not present in the rotational experiments. 

We again observed that the effective hydrodynamic 
radii obtained for the hollow and filled raspberry differ 
significantly, 3.38cr and 3.54cr, respectively. It should be 
stressed that the fact that behavior of /r£ is the same for 
both models, does not imply hydrodynamic consistency 
of the model, when we compare the value of and 
for the same model, which we will do next. 


3. The Effective Radius 


The Radius Dependence on Various Parameters 

To further assess the significance of the difference be¬ 
tween the effective hydrodynamic radii, we repeated our 
experiments for two other values of the bare friction i/o 
and several R. The results for the box-size dependence 
were in quantitative agreement. Our results for the hy¬ 
drodynamic radii are summarized in Fig. These were 
obtained, as before, by extrapolating to the bulk value 
of the mobility. The translational and rotational radii of 



Rh(F) — Rh"(H) -EH- R"(D) □ 


FIG. 8. (color online) The difference of the bulk effective 
hydrodynamic radii (translational rotational i?^) with 
respect to the imposed radius R for the filled (F), hollow 
(H), and dense shell (D) raspberry models, respectively, (a) 
The difference as a function of the bare friction coefficient 
expressed in LB units (time r and mass mo) for raspberry 
particles with a radius of i? = S.Ocr. (b) The difference as 
a function of the imposed radius R for a bare particle-fluid 
friction of = 25moT“^. 


the hollow raspberry differ substantially. This result is in 
agreement with the findings of Ollila et al. [Ml El] and it 
is in line with the theoretical predictions of Refs. |551 - HT| 
The mismatch occurs for all values of the friction coeffi¬ 
cient and radius that we examined. This discrepancy is, 
however, undesirable to simulate hard colloidal spheres, 
a purpose for which the raspberry model was initially in¬ 
troduced. m In each case, the agreement between the 
effective hydrodynamic radii of the filled raspberry par¬ 
ticles is almost perfect. 

We also performed experiments using a hollow rasp¬ 
berry with Ntot = 925 - the same total number of beads 
as in the filled raspberry - we refer to this model as the 
‘dense shell’ raspberry. This allowed us to examine the 
hypothesis that we simply obtained an increased effec- 
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live friction with the greater bead numbers used in the 
filled raspberry, leading to a better match between rota¬ 
tional and translational hydrodynamic radius. |35) Also 
note that the screening ratio for the filled and dense rasp¬ 
berry is the same, see Table |lj A discrepancy between 
and was found for the dense shell raspberry, see 
Fig.lD In fact, the deviation is slightly larger than for 
the N = 202 hollow raspberry. This can be attributed to 
an overall improvement of the coupling in the dense shell 
raspberry, which forces the translational radius towards 
the no-slip value more quickly than the rotational one. 



FIG. 9. (color online) The difference of the bulk effective 
hydrodynamic radii as a function of the number of internal 
coupling points N' for a raspberry with radius R = 3.Oct. 
The dots and connecting red curve show the translational hy¬ 
drodynamic radius RJ^ and the blue squares with connecting 
curve show the rotational hydrodynamic radius R^. 

To investigate the impact of the level of filling on the 
particle, we varied the number of internal coupling points 
N' for a raspberry with radius R — 3.Oct with N = 202 
shell-coupling points. The result is shown in Fig.[^ which 
gives the dependence of the hydrodynamic radii on the 
filling parameter N'. It is clear that the correspondence 
between and can be substantially improved by 
adding coupling points, until there is essentially no longer 
a difference, at our chosen value of N' = 722. This 
correspondence is reached at a feasible number of cou¬ 
pling points. However, it requires considerably more than 
one coupling point per lattice cell, i.e., a filling density 
greater than 10.0a“^ was found to give almost perfect 
correspondence between the two radii. 

We examined the fluid flow inside the filled and hollow 
R — 3.Oct raspberry with Ntot = 925 and 203 respectively, 
to determine the cause of the inconsistency between the 
effective hydrodynamic radii for the hollow model. Fig¬ 
ure [I^ a) shows the flow field around a hollow and filled 




r/a 


FIG. 10. (color online) Gomparison of the flow held around 
a Filed and hollow raspberry, respectively, undergoing a con¬ 
stant rotation, (a) Two dimensional plane through the center 
of the sphere with a normal that is parallel to the axis of 
rotation. The result for the Filed raspberry is shown on the 
left and for the hollow variant on the right. The color coding 
gives the magnitude of the Fuid velocity on the grid (blue 
lines). The thick black circle roughly indicates the position of 
the coupling points (at Irj = R). The dashed blue semi circle 
and half square serve as guides to the eye for the structure 
of the Fow Feld inside the raspberry, (b) Magnitude of the 
Fuid velocity Vf{r) - expressed in MD units of time, r, and 
position CT - as a function of position r along the black vertical 
divide in (a). Only the value inside of the raspberry is shown 
for the Filed (red, solid) and hollow (blue, dashed) particle. 


spherical raspberry, rotating at constant angular velocity 
about the axis pointing into the page. From the flow field 
it becomes apparent that the coupling of the raspberry to 
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the fluid has more lattice artifacts (is less smooth) for the 
hollow raspberry than for the filled one, indicating poorer 
coupling. We quantified this difference further by exam¬ 
ining the fluid velocity inside the particle, see Fig.[Toj)b). 
While the filled raspberry shows a linear increase in the 
velocity with the distance from the center (similar to the 
so-called ‘Rankine vortex state’), the hollow raspberry 
shows a clear kink in the velocity profile. This kink can 
be attributed to the diminished fluid-particle coupling 
away from the shell of MD beads. Effectively, the hol¬ 
low shell raspberry achieves its Screening ratio (a low 
Brinkman length) only close to the shell, whereas the 
filled raspberry achieves low permeability throughout. 

The Porosity of the Raspberry 



Finally, let us discuss our results in the context of the 
predictions made by theory for porous spheres. |39H41) In 
Fig. 11'a) we have plotted the theoretical prediction for 
the ratio of the hydrodynamic radii Rp/Rp as a function 
of the screening ratio kx (‘X’ is either ‘F’ for filled or ‘H’ 
for hollow). This data is based on Tableland Eqs. (27 


35). We also show the ratio R]^/R^ that we obtained 
from our raspberry simulations. It is clear that there is 
a mismatch between our results and the predictions of 
theory. That is to say, the trends predicted by theory 
are not reproduced. This can be attributed to the fact 
that the theory solves Stokes’ equation with a Stokeslet 
point-coupling. The reality of the finite grid-size LB sim¬ 
ulations is that the point-coupling only approximates the 
Stokeslet. m Correspondence is only found at a few 
lattice spacings away from the coupling point and the 
Stokeslet form can only be reproduced in the limit of 
small lattice spacings a. [12] From Fig. ©a) it becomes 
clear that for finite a, the translational hydrodynamic ra¬ 
dius is larger than that of the rotational one; the opposite 
of the theoretical prediction. [39H4T] 

This leads to the question: “Can the result of the 
theory in principle be reproduced by our simulations?” 
In order to determine this, we chose parameters which 
are unsuited to achieve our goal of obtaining hydrody¬ 
namic correspondence, but allow for a sufficiently low 
fluid-particle coupling to observe the difference in radii. 
For a hollow raspberry with radius R = 3.Oct, TV = 25 
coupling points on the surface, and a bare friction of 
Co = 5.0moT“^, the theory predicts a radius ratio of less 
than one and a strong decay of this ratio with the lattice 
’b). The curve is based on a combina- 


spacing, see Fig. 


as a“^, we used a fixed box size L = 20ct. This allowed us 
to use grid spaces of a = 0.125 ct, i.e., LB grids with 160^ 
elements, which is roughly the limit of the grid size that 
fits into a modern GPU’s memory. We therefore did not 
perform finite-size scaling. We exploited the Hasimoto 
relation |8] of Eq. (40) to fit for the effective translational 
hydrodynamic radius and the Hofman et al. |4] relation of 
Eq. (46) to fit for the effective rotational hydrodynamic 
radius. 


Fig. 

11 

(27 

35 


35). Since the computation time scales 





FIG. 11. (color online) Comparison to the results of theory for 
porous spheres. [391441) (a) The ratio of the translational and 
rotational hydrodynamic radii as a function of the screening 
ratio Kx (‘X’ is either ‘F’ for filled or ‘H’ for hollow). The 
ratio Rp /Rp from the theory |39H41 | is indicated by a red 
solid curve for a filled sphere and a blue dashed curve is for a 
hollow shell. The ratio Rf /R^ for the raspberry simulations 
are indicated using symbols with standard error bars; red dots 
show the filled raspberry results and blue open squares those 
for the hollow raspberry. The thick gray line indicates the unit 
ratio, (b) The ratios as a function of the lattice spacing a that 
were obtained for a hollow raspberry with radius R = 3. Oct, 
N — 25 coupling points on the surface, and a bare friction 
of Co = 5.0moT~^. The blue squares with error bars show 
the results of our raspberry simulations. The red solid curve 
shows the prediction of the theory for this system. The dashed 
green line shows the prediction of the theory for a slightly 
higher bare friction (Co = 8.0mor“^). 


Figure [TT](b) shows the result of our simulations. The 
error bars are sizable, but appropriate for the limited 
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box-sizes that we could study upon varying a. It is 
clear that for these parameters our data has the same 
trend as predicted by the theory. However, we found 
that a slightly higher bare friction coefficient, namely 
^0 = 8.0moT“^, yields better agreement with our sim¬ 
ulation result. This difference can be attributed to the 
fact that the theory assumes distribution of the screening 
ratio kh that is homogeneous over the shell, whereas our 
numerical results are for individual coupling points. For 
a higher number of coupling points, the decay in Rp /Rp 
is very weak for reasonable LB parameters, which makes 
it difficult to see if the theoretically predicted trends are 
matched within the error bar. An additional source of 
discrepancy is the effective shell-width of ^ 2a. [12] At 
the maximum resolution that we were able to achieve, 
there is still an effective width of 0.25 (t, whereas the 
theory assumes a dirac-delta distributed Screening ratio 
for the hollow raspberry. [39ll4T| Considering these two 
sources of error, the agreement with theory that we were 
able to achieve is quite excellent. With sufficient com¬ 
putational resources, the porosity prediction should be 
captured for a more dense distribution of shell coupling 
points and even smaller value of a, but this falls outside 
of the scope of our current investigation. 


B. Dumbbell in a Simple Cubic Crystal 

Thus far, we have concentrated on the quality of the 
raspberry approximation for convex objects, namely the 
specific case of a spherical particle. In order to assess 
the raspberry model’s ability to capture the hydrody¬ 
namic properties of a non-convex particle, we considered 
two dumbbell-shaped raspberries, as shown in Fig.|^ We 
took care to create a dumbbell raspberry model for which 
the two spheres touch, when the effective hydrodynamic 
radius of the MD beads is taken into account, see Fig. 
(left). Note that for a dumbbell-shaped particle the hy¬ 
drodynamic mobility tensor (HMT) has a diagonal form, 
with translational mobilities in the top-left 3x3 block 
(sub-matrix) and rotational ones in the lower-right 3x3 
block. There are no cross-coupling terms due to symme¬ 
try. [551170] 

Our results for the dumbbell particles are qualitatively 
similar to those shown for the spherical colloid discussed 
above. Namely, we found the box-size dependence to 
be of the form (^1 + Ai/L + Bi/, with 

X either R or T and i either T or ||, and Ai and Bi 
coefficients. However, we could not compare our re¬ 
sults to analytic calculations, since, to the best of our 
knowledge, such expressions have not been formulated 
for dumbbell-shaped particles. We therefore considered 
the extrapolated bulk mobility coefficients only. Using 
both quiescent and thermalized simulations we verified 
that the HMT has the expected form. In particular, all 
off-diagonal coefficients were orders of magnitude smaller 
than the diagonal elements and zero within the error 
bars. Moreover, we found that for both the translational 


Method 

T / T 

M|| /Mo 




d = 7a / d = 2.00 

Rasp. (H) 
Rasp. (F) 

ED 

0.78 ±0.01 
0.77 ±0.01 
0.77 ±0.01 

0.70 ±0.01 
0.69 ±0.01 
0.70 ±0.01 

0.61 ±0.01 
0.55 ±0.01 
0.55 ±0.01 

0.28 ±0.01 
0.27 ±0.01 
0.27 ±0.01 


d = 5a / d — 1.43 fj.m 


Rasp. (H) 

0.83 ±0.01 

0.75 ±0.01 

0.67 ±0.01 

0.39 ±0.01 

Rasp. (F) 

0.82 ±0.01 

0.74 ±0.01 

0.59 ±0.01 

0.36 ±0.01 

m 

0.82 ±0.01 

0.75 ±0.01 

0.60 ±0.01 

0.37 ±0.01 


TABLE II. Comparison for a dumbbell-shaped particle be¬ 
tween the results obtained using the raspberry model - both 
hollow (H) and filled (F) - and HYDROSUB/HYDRO++ [ZD 
I72| for the translational (T) and rotational (R) mobilities in 
the direction parallel (||) and perpendicular (T) to the main 
axis in bulk fluid. The mobilities are normalized by the bulk 
values for a sphere with the same radius as one of the spheres 
comprising the dumbbell. 


and rotational mobility sub-matrices, the two entries cor¬ 
responding to perpendicular motion were equal (within 
the error) and the parallel component was larger, as ex¬ 
pected. Tablejnjlists these mobility coefficients. In order 
to non-dimensionalize the results, we divided the mobility 
coefficients by the translational and rotational mobility 
of a sphere with radius R — 3.Oct (the size of one of the 
dumbbell’s lobes), respectively. 

To validate our model for the simulation of anisotropic 
non-convex particles, we compared our data with the re¬ 
sults obtained using the HYDROSUB and HYDRO++ 
program. [niiTi] These are tools used to evaluate the 
hydrodynamic properties of macromolecules and have 
been successfully utilized in comparisons to experimental 
data for solid anisotropic colloids. m We determined the 
HMT using the methods of Refs. mma for dumbbells 
consisting of two spheres with radii R = 1.0 fam at po¬ 
sitions (±1.0,0, 0) fim (touching) and (±0.714,0,0) fj,m 
(overlapping), respectively, in a fluid of viscosity 1.0-10“^ 
kgm“^s“^ and density 1.0 • 10^ kgm“^ with tempera¬ 
ture T = 293.15 K. We assumed that the particle has 
the same density as the fluid. The numerical algorithm 
is parametrized as follows: H = 26, iLmax = 1-5 ■ 10^, 
^max = 80.0 • 10“®, and A^trials = 10,000; which are 
internal commands. The number of intervals for the dis¬ 
tance distribution was set to 30. By applying the same 
numerical parameters to the case of a single sphere we 
obtained the reference data used to normalize the result, 
which in turn allows for a direct comparison to our re¬ 
sults. 

The results of this comparison are summarized in Ta¬ 
ble [H] in which we give the mobilities for the filled and 
hollow raspberry, as well as the ones determined using 
the methods of Refs. [HI El]. The agreement for the 
translational bulk mobilities is excellent in all three data 
sets. However, it is clear that for the hollow raspberry 
there is a significant difference in the rotational mobility 
ratio with respect to the result for the filled and HY- 
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DROSUB/HYDRO++ simulations. This difference lies 
well outside of the error bar of the average of the latter 
two. This confirms that our ‘filling + fitting’ procedure 
is effective for more complex (non-convex) geometries, as 
expected. 


IV. DISCUSSION 


In Section III we have demonstrated that our ‘filling + 
fitting’ formalism leads to excellent agreement between 
established theoretical and numerical results for the hy¬ 
drodynamic behavior of convex and non-convex solid par¬ 
ticles. By ‘filling -|- fitting’ one significantly improves the 
agreement between the effective hydrodynamic radii ob¬ 
tained by translational and rotational experiments, re¬ 
spectively, allowing the point-coupling LB model to de¬ 
scribe solid particles. The improvement is related to a 
reduced permeability throughout the particle - in line 
with the findings of Ref. [M] . The hollow-shell raspberry 
achieves this only locally. o EH IMl 137] In this section 
we discuss this discrepancy between the effective hydro- 
dynamic radii in more detail and place our work in the 
context of previous studies. 

The fractional difference in hydrodynamic radii of ap¬ 
proximately O.I(t/3.4(t = 0.03 for the hollow raspberry 
may seem perfectly acceptable for most applications. 
However, one should be careful, since this small frac¬ 
tion can lead to a 10% discrepancy between the expected 
translational and rotational mobility, had we assumed 
the effective hydrodynamic radius for rotational motion 
to be the same as that for translational motion. In pro¬ 
cesses involving both translation and rotation, this could 
lead to significant deviation from the desired behavior. 


Previous Studies 


A closer examination of the data presented in the orig¬ 
inal raspberry paper by Lobaskin and Diinweg m shows 
that the trends in matching to the results of Refs. HIBl- 
lllj with effective hydrodynamic radii observed in our 
work, are captured by their data points. Lobaskin and 
Diinweg erroneously assumed that the radius of the par¬ 
ticle was the same as the radius R at which they po¬ 
sitioned their MD beads. Within the numerical uncer¬ 
tainty present in their results and the computational abil¬ 
ities of the time, this extrapolation to bulk was unavoid¬ 
able. By re-examining the data points of Ref. m , we 
conclude that it is possible to fit the following bulk mo¬ 
bilities 



(48) 

= (0.90 ± 0.02) 

(49) 


This indicates that there is indeed an effective radius, 
Rl = (1.03 ± 0.02)i? and = (1.03 ± 0.0I)R, but the 
data is not of sufficient quality to assess whether there is 


a difference between the effective translational and rota¬ 
tional hydrodynamic radius in their measurements. 

Chatterji and Horbach m carried out a more thor¬ 
ough examination of the effective translational hydrody¬ 
namic radius. However, they did not provide results for 
the rotational hydrodynamic radius, they only comment 
on having carried out such experiments. Our results in 
Fig. for the value of for the hollow raspberry are in 
quantitative agreement with Ref. m- We therefore deem 
it likely that a similar discrepancy would be present in 
the data of Ref. m, especially considering our observa¬ 
tions and those of Refs. [34l 135] . 

Finally, Poblete et al. EO] did not report a differ¬ 
ence in the bulk hydrodynamic radii using their MPCD 
method for a hollow raspberry. They instead found agree¬ 
ment between the two. However, it is unclear how ac¬ 
curately Poblete et al. could extrapolate their results 
to the bulk value, as in MPCD one always works with 
thermalized and therefore noisy data. In addition, it is 
not obvious how large the effect (i?^ ^ i?^) would be 
for their high-speed-of-sound systems. Furthermore, the 
grid-shifts that are typically applied in MPCD to restore 
Galilean invariance, may substantially reduce any such 
lattice-discretization and porosity effects. 

Relation to the Work of Ollila et al. 

The inconsistency between the translational and rota¬ 
tional mobility in the raspberry model was first pointed 
out by Ollila et al. [M] [3S| Reference [33] contends that 
these inconsistencies are representative of the properties 
of the point-coupling method. Namely, that the objects 
modeled using this formalism are porous. Ollila et al. ar¬ 
gue that this porosity leads to problems when using this 
type of model to describe solid objects. In particular, 
models that fit for the radii should be considered with 
suspicion according to Ref. [33], as the fitted hydrody¬ 
namic radii may be inconsistent between various hydro- 
dynamic experiments. This assessment may seem in di¬ 
rect contradiction to our observations. However, Ollila et 
al. do not exclude the possibility of finding numerical pa¬ 
rameters for which a quality fit can be made. We have 
shown here, as well as in Ref. [7], that our ‘filling -|- fit¬ 
ting’ formalism works well to match the simulations to 
analytic results for solid particles over a wide range of pa¬ 
rameters. That is, we obtained numerically consistency 
for physically relevant hydrodynamic experiments. 

Note that the excellent agreement shown between the 
simulation results and analytic expressions for porous 
spheres in Ref. |34j is not without caveats. In partic¬ 
ular, Ollila et al. indicate that it is necessary to use 
a particle radius for the coupling points that is ‘incom¬ 
mensurate’ with the lattice to obtain the excellent corre¬ 
spondence for the translational properties of the porous 
particles without fitting. Due to the properties of the 
interpolation scheme, this incommensurability criterion 
and the subsequent choice of a particle radius that yields 
correspondence, can be treated on the same footing as a 
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fit parameter. Moreover, Ollila et al. require an effec¬ 
tive hydrodynamic radius (another fit) to obtain similar 
correspondence for the rotational properties of their par¬ 
ticles. 

We have performed our simulations with both station¬ 
ary and moving particles at positions and in directions 
both commensurate and incommensurate with the lat¬ 
tice. In all these experiments, we did not find a siz¬ 
able change in the effective radii, nor a breakdown of 
the correspondence between the two. We thus argue 
that our ‘filling -|- fitting’ method is a cleaner and more 
forthright way of proving a correspondence between a 
theoretical result and simulations. We therefore believe 
that an equally excellent correspondence between theory 
and simulations could have been achieved in Ref. |34j , by 
dropping the incommensurability criterion and fitting for 
both effective hydrodynamic radii. In addition, our ‘fill¬ 
ing -I- fitting’ method is an excellent approach to find LB 
and coupling parameters for which the behavior of solid 
objects in a Stokes’ fluid can be faithfully reproduced. 

Numerical Efficiency Considerations 

In both Refs. [SU [3S] the number of coupling points 
used to obtain correspondence between theory and sim¬ 
ulations is rather large. Such a high number of points 
is acceptable in addressing questions of a fundamental 
nature, but it may prove problematic in performing sim¬ 
ulations with many raspberry particles, as is typical for 
self-assembly and crystallization studies. m 

The algorithm may become prohibitively expensive for 
these numbers of coupling points. Lowering the over¬ 
all number of coupling points and specifically their local 
density is of particular relevance to GPU-based LB im¬ 
plementations. The force applied to the nodes of the 
LB grid by a coupling point is calculated using so-called 
‘atomicAdd’ operations. These operations can be used 
to avoid race conditions that arise from colliding mem¬ 
ory requests. However, for a large number of beads close 
to a specific LB node (high coupling-point density), the 
use of the atomicAdd operation can cause the program to 
slow down significantly. Therefore, reducing the number 
of local coupling points is of paramount importance and 
our ‘filling -|- fitting’ procedure is thus numerically favor¬ 
able to models that require a higher local coupling-point 
density. 

Porosity and the Filling Factor 

With regards to the filling procedure, we obtained rea¬ 
sonable consistency for the translational and rotational 
mobility for roughly 10 beads per LB grid volume 
within the particle, for a grid size of a = I.Oct. Depend¬ 
ing on the radius of the particle, we obtained a slight 
variation of the effective increase in hydrodynamic size 
of the particle, but for all of our R and (jo data points, 
we managed to achieve far superior agreement between 
the hydrodynamic radii for the filled raspberry than for 


the hollow variant. 

As mentioned above, in discussing the results by 
Ollila et al. [331 ES]j our results seem counterintuitive, 
when one considers theoretical predictions for porous me¬ 
dia. [55MT] Indeed, we find that for the LB parame¬ 
ters that we chose (which are physically reasonable), our 
results deviate substantially from the predicted trends. 
However, one should be careful in interpreting these re¬ 
sults. The LB method solves a discretized form of the 
Boltzmann transport equation, which only reproduces 
the result of Stokes’ equation in the appropriate limits 
(small grid spacing, etc.). This means that point-forces 
that are applied to the LB fluid via the Ahlrichs and 
Diinweg interpolated point-coupling scheme m are not 
true Stokeslets. Therefore, it is a priori not to be ex¬ 
pected that the theoretical result is reproduced. We have 
demonstrated that the results of porous sphere theory 
can in principle be recovered in our simulations, when 
the surface coverage of coupling points and the bare fric¬ 
tion coefficient are sufficiently small. Whether this result 
can be extended to greater particle numbers is difficult 
to assess with the precision we can currently achieve. 

The Short-Time Behavior 

We end with a comment on the short-time behavior of 
the raspberry particles. As originally shown in Ref. 
the Ahlrichs and Diinweg interpolated point-coupling 
scheme [13] has problems in reproducing the short-time 
properties of the (A)VACF that are expected for a solid 
no-slip particle [50] or even a porous colloid, |6^ due 
to the presence of an unphysical coupling regime. The 
correct zero-time value of C'J(O) = SksT/m is achieved, 
but there is no decay to C'[{t > 0) = 3kBT/m*, with m* 
the virtual mass, over a time scale related to the propaga¬ 
tion of sound. [50] Instead, a much lower plateau value for 
C]^{t > 0) is reached. Felderhof (50] has pointed out that 
the secondary (virtual mass) regime is not present for a 
porous colloid in the high frequency limit. It is possible 
that our findings are in agreement with this observation 
and that the decay with m* reported by Lobaskin and 
Diinweg. m is not to be expected. However, the pres¬ 
ence of the unphysical coupling regime and the difficulty 
in determining the limit in which our LB is operating 
frustrate further analysis at this time. An analysis of 
the influence of both the frequency and strength of the 
viscous coupling on the result and the correspondence to 
the predictions of Ref. jBD] is left for future study. 

It has been suggested that the modified coupling 
scheme by Mackay et al. [33] may resolve the short-time 
decay issues. However, examination of the VACFs re¬ 
ported in Ref. m reveal that the double exponential- 
type decay shown in their fluid-mass dominated result is 
not captured by the result of Ref. [SD] , as is reported in 
Ref. [37] (but not demonstrated using fitting procedures). 
This is, again, expected on the basis of the results by 
Felderhof. [50] Unfortunately, it is also not clear that the 
predictions of Ref. m are more accurately reproduced. 
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Furthermore, an analysis of the results of the oscillatory 
experiments in Ref. |34j does not yield significant addi¬ 
tional insight into the short-time quality of the Mackay et 
al. algorithm. In particular, it is unclear whether the 
Mackay et al. coupling algorithm was used in Ref. [5^ . 
In addition, the period of the oscillation is sufficiently 
long to obfuscate any short-time discrepancies that may 
be present. 

V. CONCLUSION AND OUTLOOK 

Summarizing, we have re-examined the properties of 
the hybrid LB and Langevin MD scheme for simulat¬ 
ing colloids developed by Lobaskin and Diinweg, m the 
so-called ‘raspberry’ model. We studied this model us¬ 
ing a variety of classic fluid dynamics experiments that 
predominantly focused on the long-time mobility proper¬ 
ties of these particles. We considered the hydrodynamic 
properties of spherical raspberries, as well as dumbbell¬ 
shaped raspberry particles in the low Reynolds number 
limit. Our results show that the proper solid-particle mo¬ 
bility in this limit is reproduced to a surprising degree of 
accuracy over a wide range of viscosities for both convex 
and non-convex particle shapes. 

From our combined data we can draw the following 
conclusions concerning the quality of the raspberry model 
and our ‘filling -|- fitting’ procedure to match its hydro- 
dynamic properties to that of a solid object in a Stokes’ 
fluid. 

• Using a raspberry model to approximate a parti¬ 
cle’s coupling to an LB fluid gives rise to an effective 
hydrodynamic radius. Our result is in agreement 
with the findings of Chatterji and Horbach. m 

• The short-time properties of a no-slip or perme¬ 
able colloid are not faithfully reproduced by the 
Ahlrichs and Diinweg interpolated point-coupling 
scheme, [13] as was first pointed out in Ref. [14]. 

• The traditional ‘hollow’ raspberry model - an 
empty shell of MD coupling beads that describes 
the particle’s surface - gives rise to a discrepancy 
between the translational and rotational effective 
hydrodynamic radius. This effect was first pointed 
out by Ollila et al. [SU |3S] and discussed in the 
context of porous particle dynamics. [HilMT] 

• We find that the aforementioned mismatch, when 
considered in the context of reproducing the hy¬ 
drodynamic behavior of a solid particle, can be re¬ 
duced to within an acceptable numerical tolerance 
by ‘filling’ the raspberry and ‘fitting’ for the effec¬ 
tive hydrodynamic radius. Here, we found a filling 
density of about 10 coupling points per LB grid cell 
gives satisfactory results. 

• The ‘filling -|- fitting’ procedure is not in disagree¬ 
ment with the assessment of Ollila et al. [51] that 


such a filling procedure is inherently problematic. 
Our result demonstrates that for reasonable LB pa¬ 
rameters the hydrodynamic properties of a solid 
particle can be effectively matched and to within 
a far higher tolerance than is possible for the hol¬ 
low variant. 

• Moreover, the formalism is not in contradic¬ 
tion with the theoretical results for porous parti¬ 
cles. [5^MT] The fact that we can obtain excellent 
correspondence between the hydrodynamic radii 
for a filled raspberry and we find inferior agreement 
for the hollow variant, should be considered in the 
context of the level of discretization that is used. 
Hence, a mismatch between the theoretical predic¬ 
tions and the numerical result is not unexpected; 
in this case it can be exploited. 

• The ‘filling -|- fitting’ procedure can be used to 
improve the raspberry model’s ability to simulate 
both convex and non-convex solid particles. We 
verified this for the specific case of a dumbbell¬ 
shaped particle, and our procedures may be safely 
extrapolated to more complicated shapes. 

• The result of Ollila et al. [S^ suggest that a regime 
can be found for which the hydrodynamic hull is 
sufficiently shrunk that it matches with the im¬ 
posed position of the coupling points. However, 
a prohibitive number of coupling points may be re¬ 
quired to achieve this condition. This is especially 
problematic for GPU-based algorithms. Our ‘filling 
-I- fitting’ procedure allows us to use a substantially 
reduced number of coupling points and still obtain 
excellent numerical agreement. 

• The force and velocity experiments traditionally 
performed to determine the translational mobility 
in a cubic geometry with periodic boundary condi¬ 
tions are problematic for small boxes compared to 
the particle size. The back force/velocity density 
that must be applied to the fluid to maintain zero 
center of mass velocity, leads to difficulties in in¬ 
terpreting the mobility data that is obtained from 
these experiments. A possible solution to this prob¬ 
lem is to identify the node locations at which the 
particle is found and to only apply the properly- 
rescaled counter force elsewhere. 

From the above, it becomes clear that the raspberry 
model is an excellent way to approximate the long-time 
regime of the fluid-particle coupling for a solid object in 
an LB fluid. However, there remains several open prob¬ 
lems to be addressed in future studies. We have shown 
that the short-time behavior of the raspberry model (for 
the LB parameters used in this manuscript) is quite dif¬ 
ferent from the low Reynolds number solution to Stokes’ 
equations. [50] |60] This raises the question of how ac¬ 
curately the short-time regime of colloid dynamics can 
be captured using the raspberry or any point-coupling 
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model. A faithful reproduction of such short-time pro¬ 
cesses would be relevant for, e.g., nucleation and crystal¬ 
lization. m Despite this concern, our analysis stresses 
the power of our ‘filling -|- fitting’ method as a means 
to find parameters for which the translational as well as 
rotational hydrodynamic properties of the original rasp¬ 
berry model can be sufficiently enhanced to give rise to 
Stokesian fluid - solid particle coupling. 
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